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Abstract — A  three-point,  sixth-order,  staggered,  combined  compact  difference  (SCCD)  scheme 
is  proposed  for  numerical  models.  The  SCCD  scheme  is  a  generalization  of  the  previously  proposed 
combined  compact  difference  (CCD)  scheme,  and  has  major  improved  features  such  as  error  and 
computational  (CPU)  time  reduction  especially  for  odd-order  difference  equations  with  odd-number 
boundary  conditions.  For  nonperiodic  boundaries,  an  additional  sixth-  or  fifth-order  boundary  condi¬ 
tion  is  proposed.  The  stability  of  the  SCCD  scheme  is  studied  using  the  eigenvalue  analysis.  ©  2000 
Elsevier  Science  Ltd.  All  rights  reserved. 
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1.  INTRODUCTION 

Following  the  trend  toward  highly  accurate  numerical  schemes  of  partial  differential  equations 
(PDE)  led  by  many  authors  (e.g.,  Adam  [1];  Lele  [2];  Carpenter  et  al.  [3,4];  Chu  and  Fan  [5]; 
Hirsh  [6];  Rubin  and  Khosla  [7];  Navon  and  Riphagen  [8],  and  Chu  and  Fan  [9,10])  proposed 
a  three-point  sixth-order  uniform  and  nonuniform  combined  compact  difference  (CCD)  schemes 
to  increase  accuracy.  This  scheme  follows  earlier  work  on  use  of  second  derivatives  in  compact 
difference  (such  as  [7]) 

l 

Y.  («*/»+*:  +  hfl+k  +  cfc/i+fc)  =  (!) 

fc=-i 

which  is  referred  to  as  the  Hermite  formula.  Here,  /  is  a  dependent  variable.  Lele  [2]  and 
Carpenter  et  al.  [3,4]  performed  detailed  studies  of  several  compact  schemes  in  terms  of  accuracy 
and  efficiency,  both  for  first-  and  second-order  derivatives  f  and  f".  Moreover,  Carpenter  et 
al.  [3]  address  the  issue  of  stability  of  boundary  stencils  from  the  theoretical  point  of  view. 

Many  ocean  and  atmospheric  numerical  models  are  designed  on  staggered  grids,  this  is  due 
to  the  fact  that  the  staggered  grids  give  more  accurate  wave  propagation.  Applying  the  CCD 
scheme  to  the  staggered  grids,  we  obtain  a  three-point  sixth-order  staggered  combined  compact 
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difference  (SCCD)  scheme  with  sixth-order  or  fifth-order  accuracy  at  both  interior  and  boundary 
points. 

2.  SCCD  SCHEME 

2.1.  General  Form 

Let  /  ( x )  be  defined  on  the  interval,  0 <  x  <  L.  We  discretize  the  interval  into  0  =  x0  <  x\  < 
•  •  •  <  xn  =  L  for  nonperiodic  boundaries  (Figure  la)  and  into  0  <  xq  <  xi  <  ■  ■  •  <  xn  <  L 
for  periodic  boundaries  (Figure  lb)  with  a  grid  spacing  h  =  Xi  -  x*_i  =  L/N,  (i  —  1,2, ... ,  N). 
Here,  Xq,  x\,  . . .  ,£jv  are  unstaggered  points. 


(b)  Periodic  boundaries. 

Figure  1.  Illustration  of  the  SCCD  grid  structure. 


Use  staggered  grids,  0  =  xso  —  xo  <  xsi  <  X\  <  •••  <  xsn  <  %n  =  %sn+ l  =  L  for 
nonperiodic  boundaries  (Figure  la),  and  0  =  xsi  <  xi  <  xs2  <  •••  <  xsn  <  xn  <  L  for  periodic 
boundaries  (Figure  lb)  with  xsi  =  -  h/ 2,  i  =  1,2,  3,  ...,1V.  Here,  xsi,  %S2,  •  ■  ■  ,%sn  are 

staggered  grid  points.  Notice  that  staggered  and  unstaggered  points  are  collocated  at  nonperiodic 
boundaries  (xso  =  %o  and  xsn+ l  =  x^). 

The  dependent  variable  f(x)  and  its  second  derivative  f"(x)  are  given  at  the  unstaggered 
points.  The  first  derivative  f'(x)  is  given  at  the  staggered  points.  By  application  of  the  Hermite 
formula  (1)  to  the  staggered  grids,  we  obtain  a  general  form  of  the  SCCD  scheme 

an/si-i  +  fsi  +  ai2/i+i  +  Pnfi-i  +  Puf”  =  7i(/i  ~  /t-i),  (2) 

a2lfsi  +  Ol22fsi+l  +  P2lfi-!  +  fi  +  p22fi+l  =  72l/i-l  —  (721  +  722  )/i  +  722/t+l-  (3) 

In  computing  first  and  second  derivatives,  systems  (2)  and  (3)  will  be  solved  for  unknown  variables 
/'  and  f"  as  /  is  given  at  the  regular  grids  fi- 1,  ft,  and  fi+\.  The  coefficients  in  (2)  and  (3)  are 
determined  by  the  Taylor  series  expansion. 
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2.2.  Expansion  on  Staggered  Grid  Points 
2.2.1.  General  expressions 

Expansion  of  /,  /',  and  f"  into  the  Taylor  series  at  the  staggered  grids  x$i  and  elimination  of 
h ?  and  h4  terms  leads  to  (Figure  2a) 

f's,  -  it  (&-,  +  f'sW )  +  ig  (/"  -  //l.)  =  ^  (/.  -  (4) 

with  truncation  error  —  l/(7!)(457/2032)/^)/i6.  Using  a  similar  procedure  expansion,  we  also  can 
obtain 

fSi  =  2  (/<  +  fi-i)  +  32^  (fsi+i  ~  f'si-i)  ~  (f?  +  f>ij)  ,  (5) 

with  truncation  error  -l/(6!)(137/27)/gp/i6. 
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(b)  Standard  grid  points. 

Figure  2.  Illustration  of  the  SCCD  grid  structure  at  interior. 

In  some  cases,  if  the  third-order  derivative  is  required,  a  similar  formulation  will  be 

fHi  —  __2AL_  ( t'_  _l  t'  \  i  ^37  ( tii  t"  \  i  240  /  j.  j,  N 

JSi  127/)2  v7St-l  +  JSi+l)  '  1971,  \Ji  Ji-l)  '  -,97.3  (./«  Ji-l)  ■ 


bt  127h2  •,Si+i;  1 27h  yJi  ~  Ji~l>  127 u<  ~  /i-u  •  W 

Expressions  (4)-(6)  are  valid  from  x$i  to  x$n  for  the  periodic  boundaries,  and  from  Xs2  to 
xSN-i  for  the  nonperiodic  boundaries.  Thus,  each  of  (4)-(6)  contains  N  algebraic  equations  for 
the  periodic  boundaries  and  ( N  —  2)  algebraic  equations  for  the  nonperiodic  boundaries. 

2.2.2.  Special  expressions  for  nonperiodic  boundaries 

Expansion  of  /,  /',  and  f"  into  the  Taylor  series  at  grids  x$i  from  neighboring  points  and 
elimination  of  h2  and  h4  terms  lead  to  (Figure  3) 

JLf'  +  Vf'  _A +  57  f  .32/!  7 

135/o+25/51  675 7s2  45/o  +  25 50//°  +  25h  50//2  ^ 


for  the  left  boundary  with  the  truncation  error  as  l/(7!)(2/5 )fgh6.  Similarly,  we  have 

+Ef’  _224#'  h,n  Oh  _  57  32/jv_i  7 

135/n+  25/s*  675-'s^-1  +  45/jv  25 ‘'*-1  50//*  25h  +  50//*~2  ^ 

for  the  right  boundary  with  the  same  order  of  the  truncation  error. 
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Figure  3.  Illustration  of  the  SCCD  grid  structure  at  the  left  boundary. 

Besides,  the  sixth-order  interpolation  near  the  boundary  staggered  grid  points  are 
4455 


fs  l  = 


3584 


f  _V_f  ,  3  45  9 

112  ^  bW2  +  256-'0  +  7lsi  +  1792-'0  +  896 


and 


_  4455  27  1  ’  81  h  3 h  ,  45 h2  „ 

*SN  3584 *N  \WN~l  512^-2  256 ^sw+1  7  ^  +  1792^  +  896^_1' 


The  third-order  derivative  at  the  neighboring  grid  points  will  be 

80 
9/i2 


15  (f  M  ,  25  80  _/£  6/£ 

JSl  4^3  (72  M+wh3fo  Qh2JS2  6h  +  h  . 


15  ..  f  ,  .  25  80  /ft  6/^_j 

JSAf  -  WJV-2  -  /ATj  +  -  g^2 /sat— 1  + 


4/j.3  V^V-"  •/iV/  ’  18/l2 

2.3.  Expansion  on  Unstaggered  Grid  Points 


(9) 

(10) 

(11) 

(12) 


Expansion  of  /,  /',  and  f"  into  the  Taylor  series  at  the  unstaggered  grids  x*  and  elimination 
of  h'2  and  h4  terms  leads  to  (Figure  2b) 

(fsi+ 1  -  fsi )  (/i-i  +  /,+i)  =  (fi- 1  -  2/i  +  /m)  (13) 


94 


47/i2 


with  truncation  error  — l/(7!)(l/4)/^8^  (xj)/i6.  Expression  (13)  is  valid  from  x\  to  xjy  for  the 
periodic  boundaries,  and  from  xi  to  x^r-i  for  the  nonperiodic  boundaries.  Thus,  equation  (13) 
contains  N  algebraic  equations  for  the  periodic  boundaries  and  (N  —  1)  algebraic  equations  for 
the  nonperiodic  boundaries. 

2.4.  Periodic  Boundaries 
For  periodic  boundaries 

/o  = /at>  /at+i  = /1,  fso  —  fsN’  fsN+i  ~  /si»  /o  = /ft>  and  f'bl+l  =  /”)  (14) 

the  SCCD  scheme  (4)  and  (13)  automatically  provide  the  sixth-order  accuracy  for  whole  region 
(interior  and  boundary  points). 

3.  ALGORITHM  FOR  NONPERIODIC  BOUNDARIES 

Treatment  for  nonperiodic  boundaries  is  an  important  task  for  the  high-order  schemes.  Follow¬ 
ing  Carpenter  et  al.  [11],  Chu  and  Fan  [9,10]  proposed  a  boundary  treatment  for  the  three-point 
sixth-order  unstaggered  CCD  schemes.  Here,  we  propose  a  fifth-order  boundary  algorithm  for 
SCCD  scheme. 
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3.1.  Finite  Difference  Equation  (FDE) 

Consider  a  one-dimensional  Euler  equation 

df  df 

y  +  =  d,  0  <x<L,  0  <t  (15) 

with  the  initial  and  boundary  conditions 

f{0,x)  =  g(x),  f(t,0)=p(t).  (16) 

Use  of  the  SCCD  scheme  (5),  (9),  and  (10)  leads  to  the  discretization 

^L+cf'Si  =  d,  *  =  0, 1, . . .  ,iV  +  1,  (17) 

which  contains  ( N  +  2)  algebraic  equations. 

There  are  (3 N  +  4)  unknown  variables  in  this  system:  fi  (i  =  0,  1, . . . ,  N),  fSi  (i  —  0,  1,2,..., 
N  4-  1),  and  f"(i  =  0,  1, . . . ,  N),  however,  only  (3 N  4-  2)  algebraic  equations  are  available  in  (4), 
(7),  (8),  (13),  and  (17).  We  need  two  more  algebraic  equations  to  close  the  system.  The  two 
additional  boundary  conditions  for  the  SCCD  scheme  are  (Figure  3) 

/o  +  y/si-y/o  +^/"  =  ^(/i-/o)  +  ^(/2-/o),  (18) 

f n  +  —f'sN  +  y/w - y/w-i  =  ~y  (/w- 1  ~  /at)  ~  ^  (/at-2  -  /at)  (19) 

with  a  fifth-order  truncation  error,  —  l/(180)/i^/i5,  for  the  left  and  right  boundaries,  respectively. 
The  coefficient  matrix  is  a  diagonal  system  with  eight  bands  when  nonperiodic  (Figure  4)  and 
with  seven  bands  when  periodic  (Figure  5).  Thus,  the  memory  requirement  is  8(3 TV  +  4)  when 
nonperiodic  and  7(3 N  +  4)  when  periodic. 


from  eq. 


Figure  4.  Coefficient  matrix  of  the  SCCD  scheme  for  finite  difference  calculation 
with  nonperiodic  boundaries. 
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from  eq. 


Figure  5.  Coefficient  matrix  of  the  SCCD  scheme  for  finite  difference  calculation 
with  periodic  boundaries. 


3.2.  Finite  Difference  Calculation 

The  SCCD  can  also  be  used  to  calculate  the  high-order  finite  difference  as  the  explicit  scheme. 
There  are  (21V+3)  unknown  variables:  f'Si  (i  =  0, 1, 2, . . . ,  N+ 1)  and  f"  (i  =  0, 1, 2, . . . ,  N)  with 
only  (2 N  +  1)  equations  (4),  (7),  (8),  (13),  (18),  and  (19).  Two  additional  boundary  conditions 
are  needed.  We  propose 

K  -  f  ■ +  m;* + ^  - 1 (/l  -  M  -  § 1 [h  - « + sp  <*  -  M  (20) 

for  the  left  boundary  and 


f/t  _  56  ri  288  ,  267 

-  ~tN- 1  -  J^fsN+l  ~  JSN  ~  ^2 

21  7 

~Sh?  ^N~2  ~  ^  +  45/?  U N~ 3  _  f”) 


for  the  right  boundary  with  a  fifth-order  truncation  error  —  109/(16800)/i^/i6.  The  finite  differ¬ 
ence  calculation  becomes  to  solve  2 N  +  3  equations  consisting  of  (4),  (7),  (8),  (13),  (18)— (21)  to 
get  2N  +  3  unknown  variables;  f'Si  (i  =  0,  1, . . . ,  N  +  1)  and  /"  (i  =  0, 1, . . . ,  N).  The  coeffi¬ 
cient  matrix  is  a  pentagonal  system  for  both  nonperiodic  (Figure  6)  and  periodic  boundaries,  see 
Figure  7.  Thus,  the  memory  requirement  is  5(2 N  +  3). 


1 


4.  STABILITY  ANALYSIS 

Gustafsson  et  al.  [12]  developed  stability  analysis  technique  based  on  normal  modal  analysis. 
Their  work  (generally  referred  to  as  GKS  theory)  established  conditions  the  inner  and  boundary 
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from  eq. 


Figure  6.  Coefficient  matrix  of  the  SCCD  scheme  for  solving  finite  difference  equation 
with  nonperiodic  boundaries. 


from  eq. 


Figure  7.  Coefficient  matrix  of  the  SCCD  scheme  for  solving  finite  difference  equation 
with  periodic  boundaries. 


schemes  must  satisfy  to  ensure  stability.  Later  on,  Lele  [2]  and  Carpenter  et  al.  [3]  improved  the 
GKS  technique  and  established  a  standard  yardstick  for  acceptance  of  a  new  compact  high-order 
scheme. 

We  take  the  Euler  equation  (15)  with  initial  and  boundary  conditions  (16)  as  an  example  for 
numerical  stability  analysis  on  the  SCCD  scheme  using  the  eigenvalue  analysis  [2,3].  Define  the 
error  function  by 

®rr  —  f  /> 

where  /  and  /  are  the  exact  and  numerical  solutions  of  the  Euler  equation  (15)  with  the  initial 
and  boundary  conditions  (16).  The  stability  of  the  numerical  solution  for  (15)  is  the  same  as  for 


330 


P.  C.  Chu  and  C.  Fan 


the  error  equation 

=0,  0  <  x  <  L,  t>  0  (22) 

at  ox 

with  the  boundary  condition 

err(t,  0)  =  0.  (23) 

The  domain  is  discretized  into  N  equal  intervals  Ax  with  iV+1  grid  points  x*  (i  —  0, 1, 2, . . . ,  N) 
and  the  SCCD  scheme  is  posted  at  the  N  +  2  stagged  grids  (si,  i  =  0,  1, ...  ,N  +  1)  (Figure  1). 
We  reorganize  (7),  (4),  and  (8)  into 


rewrite  (18),  (13),  and  (19)  into 
and  rearrange  (9),  (5),  and  (10)  into 


Axe'  +  Bxi"  =  Cxe, 

(24) 

A2e!  +  B2e”  =  C2e, 

(25) 

e9  =  A3e'  +  B3e", 

(26) 

where 


e  = 

e'rr{s0)Ax 
err  (#1) 
err  ( X2 ) 

,  e'  = 

e'rr(8l) 
e'rr  (s2) 

err  (xn) 

.e'rr  (sn+i)Ax_ 

e'rr  {SN-1 
.  Kt(sn) 

["  err  (xo)  1 

&rr  (^l) 

e"r(xl) 

Crr  (^2) 

e"  = 

-err  (xn). 

Ax2 

es  = 

CO 

u 

CD 

Ax  (AT,  N) ,  Bx(N  +  l,N),  Cx(N  +  2,N),  A2(N,N  +  1), 
B2{N  +  l,N  +  l),  C2(N  +  2,N  +  1),  A3(N,N),  B3(N  +  1,N ) 

are  coefficient  matrices.  Elimination  of  C\  and  C2  from  (24)-(26)  leads  to 


e!  —  Be,  e 9  =  Ae, 


(27) 


where 

B={A1-BxB^1A2)~1  (Cx-BxB^Ci),  A  =  A3B  +  B3  (B^ C2  -  B^1  A2B)  .  (28) 


Substitution  of  (27)  into  (22)  leads  to 


(29) 


Using  the  normal  modes 


e(x,t)  =  exp  (wt)e(x). 
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Figure  8.  Eigenvalues  cj'  of  equation  (30)  for  various  horizontal  resolutions. 


Equation  (29)  becomes  an  eigenvalue  problem 

ujAx 


Ae  =  Be, 


(30) 


where  e(x)  is  an  (N  +  2)-dimensional  vector.  The  algebraic  equation  (30)  has  (N  +  2)  eigenvalues 

j  =  0,1, . . .  ,N  +  1. 


,  u)jAx 

u'i  =  — - , 

3  c 


If  the  real  part  of  all  eigenvalues 

Re  (a >j)  <0,  j  =  0, 1, . . . ,  N  +  1, 


(31) 

(32) 


the  numerical  schemes  are  stable. 

We  use  the  Matlab  [13]  built-in  function  ‘eig’  to  compute  the  eigenvalues  of  (30).  Figure  8 
shows  the  eigenvalues  {w'  }  in  the  complex  plane  with  different  resolution  N.  All  Re  (w' )  are 
negative,  which  indicates  the  stability  of  the  SCCD  scheme. 
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(33) 


5.  EXAMPLES 

5.1.  Boundary  Layer  Flow  Near  a  Flat  Plate 

5.1.1.  Blasius  equation 

Laminar  flow  near  a  flat  plate  leads  to  the  Blasius  equation 

?"  +  f  f"  =  0,  0  <  x  <  10 

with  the  boundary  conditions 

/( 0)=0,  /7(0)  =  0,  //(10)  =  1.  (34) 

This  equation  does  not  have  the  exact  solution,  but  its  second  derivative  f"  has  the  exact  value 
at  the  left  boundary  [14], 

f"  (0)  =  0.469600.  (35) 

The  FDE  for  (33)  is  written  by 

/s"  +  fsj's,  =  0,  *  =  1,2, ...  ,7V.  (36) 

We  use  this  dynamical  system  to  test  the  accuracy  and  CPU  time  of  the  SCCD  scheme. 

5.1.2.  SCCD  scheme 

The  SCCD  FDE  requires  f"  on  the  staggered  grids.  Thus,  we  need  to  interpolate  /  and  f" 
into  the  staggered  grids  for  (33), 


13 


11 


277  c,  4 


t"  _  u  t"  i  t"  f  f 

/si  -  ^/o  +  jg/l  -  ^/<>  -  3^/52  - 


231 


39 


64  h2*0  +  h2^  +  64/i2 


h  > 


fsi  =  2  (/"  +  f"-i)  -  ^  (fsi+i  ~  fsi- 1) .  i  =  2, 3, . . . ,  TV  -  1, 


(37) 

(38) 


13 


11 


277  „ 


4 


231 


39 


f*N  __32/^  +  l6/^-1  +  96hf'SN+1  +  W3"-1  ~  m?fN  +  h?fN~l  +  64h2/w-2,  (39) 
Substitution  of  (5),  (6),  (9)— (12),  (37)-(39)  into  (36),  and  together  with  (4),  (7),  (8),  (13),  (18), 
and  (19)  lead  to  (37V+4)  algebraic  equations  for  (37V+4)  variables:  /,  (i  =  0, 1, 2, . . . ,  TV),  fSi  ( i  = 
0, 1, ....  TV  4- 1),  and  //'  (i  =  0, 1, . . . ,  TV). 

5.1.3.  Second-order  scheme 
The  second-order  central  staggered  scheme  is  given  by 

fi  +  fi-l 


fsi  = 

fi  = 


fi-l  ~  %fi  +  /t+1 

h? 


*  —  1,2, ...  ,7V, 
i  =  1, 2, ....  TV  -  1. 


The  Blasius  equation  (33)  is  discretized  by 

£H  ft!  fff  I  r/t 

Ji  J  i—  1  ,  r  i/i— 1  n 

- h - +fs‘ - 2 - °' 

and  the  boundary  conditions  (34)  is  changed  into 


*  =  1,2,  ...,7V 


fo  =  0,  ft  = 


«fi-fo)/h-ftt  fi-fo 


h/2 


h? 


cn  _  (fN  ~  ( fN  -  fN-l)fh)  (1  -  ( fN  -  fN-l)/h) 

Jn  h/2  h/2 


(40) 

(41) 

(42) 

(43) 

(44) 


The  Blasius  equation  (33)  with  the  boundary  condition  (34)  becomes  a  set  of  nonlinear  algebraic 
equations  (40)-(44)  for  /,,  /"  (i  =  0, 1, 2, . . . ,  TV),  and  fSi  (i  =  0, 1, 2, . . . ,  TV  +  1). 
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5.1.4.  Iteration 

The  iterative  method  is  used  to  solve  the  nonlinear  algebraic  equation  (36).  For  each  step  of 
iteration  m,  the  variables  fs{  (i  —  0, 1, 2, . . . ,  N  +  1)  are  taken  the  values  of  the  previous  step, 
Thus,  the  two  sets  of  nonlinear  algebraic  equations  become  linearized  and  easy  to  solve. 
The  iteration  starts  from  a  linear  field, 


fi0)  =  %i . 


i  =  1,2, ...  ,iV, 


which  provides  the  first  guess  values  for  fs{ 


r(0)  _  "t"  %i—  1 

J Si  ~  o  5 


<  =  1,2 . iV. 


(45) 


The  solutions  of  the  two  linearized  systems  lead  to  a  set  of  new  values  /|V ,  which  are  used  for 
the  next  iteration.  The  iteration  stops  at  the  step  M  as 


|/"<">(0)  -  /"(M-1>(0)  <  10-9. 


(46) 


5.1.5.  SCCD  verification 

We  solve  the  Blasius  equation  (33)  numerically  with  both  SCCD  and  second-order  schemes 
under  various  horizontal  resolutions  and  record  the  CPU  time  (a  SGI  origin-100  was  used).  Since 
f"  of  the  Blasius  equation  has  the  exact  value  at  the  left  boundary  (35),  the  accuracy  evaluation 
between  SCCD  and  second-order  schemes  is  pursued  using  /"( 0).  The  difference  between  the 
numerical  and  exact  values  of  /"( 0)  is  regarded  as  the  computational  error. 

Error  versus  grid  cell  number  diagram  (Figure  9a)  shows  that  for  the  same  grid  cell  number,  the 
error  of  the  SCCD  scheme  is  more  than  three  orders  of  magnitude  smaller  than  the  second-order 
scheme.  For  example,  the  errors  of  the  SCCD  and  the  second-order  schemes  for  the  cell  number  50 
are  1.28xl0-4  and  1.99xl0_1,  an  error  reduction  by  1555  times  using  the  SCCD  scheme.  For 
the  same  accuracy,  the  SCCD  requires  much  less  (1-2  order  of  magnitudes)  resolutions.  For 
example,  the  grid  cell  number  should  be  10,000  (second-order  scheme),  or  85  (SCCD)  if  the  error 
is  0.3xl0-7. 


Figure  9.  Error  and  CPU  time  comparison  between  SCCD  (solid  curve)  and  second- 
order  (dashed  curve)  schemes. 
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ICf6  10‘4  10'2 

Error 

(b)  CPU  time  versus  error. 

Figure  9.  (cont.) 

Usually,  the  CPU  time  increases  with  the  accuracy.  CPU  time  versus  error  diagram  (Figure  9b) 
shows  a  more  rapid  increase  using  SCCD  than  the  second-order  scheme,  and  the  existence  of  a 
critical  error  (5.01  x  10-5),  where  both  SCCD  and  second-order  schemes  consume  the  same  CPU 
time  (near  0.0088  s).  For  the  error  less  (more)  than  this  critical  value,  the  SCCD  scheme  consumes 
more  (less)  CPU  time. 

5.2.  cr-Coordinate  Ocean  Model 

5.2.1.  cr-coordinate 

In  coastal  ocean  prediction  models  (most  staggered),  the  effects  of  bottom  topography  must  be 
taken  into  account.  This  can  be  done  by  using  a  terrain-following  cr-coordinate  system,  where  the 
water  column  is  divided  into  the  same  number  of  grid  cells  independence  of  depth.  Let  (x„,  y,,  2) 
denote  Cartesian  coordinates  and  (x,  y,  a)  sigma  coordinates.  In  most  sigma  coordinate  ocean 
models,  the  relationship  between  the  two  coordinate  systems  are 

x  =  x.,  »  =  ¥.,  *=1  +  2^^,  (47) 

where  z  and  a  increase  vertically  upward  such  that  z  —  0,  <r  =  1  at  the  surface  and  z  =  -H  and 
a  =  —  1  at  the  bottom.  H  =  H(x,  y)  is  the  bottom  topography. 

5.2.2.  Major  problem  of  cr-coordinate 

A  problem  has  long  been  recognized  in  computing  the  horizontal  pressure  gradient  in  the 
cr-coordinate  system  (e.g.,  [5,15-17]):  the  horizontal  pressure  gradient  becomes  a  difference  be¬ 
tween  two  terms,  which  leads  to  a  large  truncation  error  at  a  steep  topography.  How  to  reduce  the 
horizontal  pressure  gradient  error  is  a  key  issue  when  using  cr-coordinate  ocean  models,  especially 
of  using  primitive  equation  models  for  coastal  regions. 

We  restrict  our  attention  to  two  dimensions  (x,  cr)  for  simplicity.  Horizontal  gradients  in  the 
2-  and  c7-coordinates  are  related  by 

dp  _  dp  .  .  dH  1  dp 

dx*  dx  a  dx  H  da 


(48) 
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Most  ocean  models  are  hydrostatically  balanced,  i.e., 


(49) 


where  g  is  the  acceleration  due  to  gravity,  p  is  the  density  of  the  fluid, 
into  (48)  and  use  of 


d  2d 


dz  H  da 


Substitution  of  (49) 


yield 


dp 

dx* 


dp 

dx 


1,.  v  dH 

j(l -»)»■£•. 


(50) 


which  is  used  for  computing  the  horizontal  pressure  gradient  in  the  cr-coordinate  ocean  models. 
Since  horizontal  gradients  -§*r,  ,  and  are  evaluated  at  the  staggered  points,  we  have  to  use 

(5),  (9),  and  (10)  to  interpolate  p  at  the  staggered  points. 


5.2.3.  Seamount  test  case 

Chu  and  Fan  [5]  used  an  ordinary  seven-point  sixth-order  difference  scheme  to  reduce  the 
horizontal  pressure  gradient  error.  Since  most  ocean  numerical  models  use  three-point  staggered 
schemes,  it  is  hard  to  apply  the  ordinary  sixth-order  schemes  (more  than  three  points)  to  ocean 
models.  For  staggered  grid  ocean  models,  the  SCCD  scheme  overcomes  the  grid  mismatch.  It  is 
easy  to  apply  SCCD  to  ocean  models. 

We  use  a  seamount  test  case  to  illustrate  the  benefit  of  using  SCCD.  Suppose  a  seamount 
located  inside  a  periodic  /-plane  (/o  =  10-4s-1)  channel  with  two  solid,  free-slip  boundaries  along 
constant  y.  The  domain  is  a  periodic  channel,  320  km  long  and  320  km  wide.  The  channel  has  a 
far-field  depth  ??max  and  in  the  center  includes  an  isolated  Gaussian-shape  seamount  (Figure  10) 
with  a  width  W  and  an  amplitude  rjs, 

v(x,  y)  =  r?max  ~  Vs  exp  - - —  - —  ,  (51) 

where  (xo>2/o)  are  the  longitude  and  latitude  of  the  seamount  center  [5]. 


Y  (km)  _15°  -150  X  (km) 


Figure  10.  Seamount  geometry. 
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The  fluid  is  exponentially  stratified  and  initially  at  rest.  The  initial  density  field  has  the  form, 


Pi  =  p(z)  +  pexp 


(52) 


where  z  is  the  vertical  coordinate,  and  Hp  =  1000  m,  and 

P{z)  =  28  —  2  exp  (53) 

is  a  reference  density  field.  Here  a  constant  density,  1000  kg  m-3  ,  has  been  subtracted  for  the 
error  reduction.  For  a  standard  test  case  [5],  the  parameters  are  given  by 

W  =  40km,  ?7max  =  5000m,  t?3  =  4500  m,  p  =  0.2kgm-3,  Ah  =  1010  m4  s-1.  (54) 

Unforced  flow  over  seamount  in  the  presence  of  resting,  level  isopycnals  is  an  ideal  test  case 
for  the  assessment  of  pressure  gradient  errors  in  simulating  stratified  flow  over  topography.  The 
flow  is  assumed  to  be  reentrant  (periodic)  in  the  along  channel  coordinate  (i.e.,  x-axis).  This 
seamount  test  case  has  a  known  solution,  i.e.,  zero  velocity,  after  integrating  from  rest  with  no 
external  forcing  (atmospheric  or  lateral  boundary)  if  there  is  no  computational  error.  Another 
nonzero  velocity  generated  is  due  to  the  computational  error.  We  use  this  seamount  case  of  the 
semispectral  primitive  equation  model  (SPEM)  version  3.9  to  test  the  new  difference  scheme.  The 
reader  is  referred  to  the  original  reference  [18]  and  the  SPEM  3.9  User’s  Manual  [19]  for  detailed 
information.  In  the  horizontal  directions,  the  model  uses  the  staggered  C-grid,  second-order  finite 
difference  discretization  except  for  the  horizontal  pressure  gradient,  which  the  user  has  choice  of 
either  second-order  or  fourth-order  difference  discretization  [20].  In  the  vertical  direction,  the 
model  uses  a  boundary  fitted  cr-coordinate  system.  The  discretization  is  by  spectral  collocation 
using  Chebyshev  polynomials.  Our  model  configuration  is  similar  to  that  of  Beckmann  and 
Haidvogel  [21],  The  time  step  and  grid  size  used  here  are 


At  =  675  s,  Ax  =  Ay  =  5  km. 


(a)  Peak  error  velocity. 

Figure  11.  Peak  error  velocity  propagation  in  20  days  for  the  SCCD,  second-order, 
fourth-order,  sixth-order  ordinary,  and  sixth-order  compact  schemes  (all  staggered). 
The  formulas  for  the  compact  schemes  compared  are  listed  in  the  Appendix. 
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Time  (days) 

(b)  Relative  peak  error  velocity  to  SCCD  scheme. 

Figure  11.  (cont.) 

Owing  to  a  very  large  number  of  calculations  performed,  we  discuss  the  results  exclusively  in 
terms  of  the  maximum  absolute  value  the  spurious  velocity  (called  peak  error  velocity)  generated 
by  the  pressure  gradient  errors.  Figure  11a  shows  the  time  evolution  of  the  peak  error  velocity 
for  the  first  20  days  of  integration  with  six  different  schemes.  The  formulas  of  these  schemes 
are  given  in  the  Appendix.  The  peak  error  velocity  fluctuates  rapidly  during  the  first  few  days’ 
integration.  After  five  days  of  integration,  the  peak  error  velocity  show  the  decaying  inertial 
oscillation  superimposed  into  mean  values:  0.4 cm/s  (ordinary  second-order  central,  three-point), 
0.048  cm/s  (ordinary  fourth-order  with  five-point),  0.04  cm/s  (compact  fourth-order  with  four- 
point),  0.01  cm/s  (ordinary  sixth-order  with  seven-point),  0.004 cm/s  (compact  sixth-order  with 
five-point),  and  0.0013 cm/s  (SCCD  with  three-point).  The  steady  approach  of  the  peak  error 
velocities  to  these  values  for  the  six  schemes  indicates  the  stability  of  the  computation.  Also,  we 
computed  the  error  ratio  between  the  existing  schemes  versus  the  SCCD  scheme  (Figure  lib). 
The  errors  for  the  second-order  (central),  ordinary  fourth-order,  compact  fourth-order,  ordinary 
sixth-order,  and  compact  sixth-order  schemes  are  roughly  300,  36.6,  30,  8,  3  times  of  the  that  of 
the  three-point  sixth-order  SCCD  scheme. 

6.  CONCLUSIONS 

(1)  From  this  study,  it  can  be  stated  that  the  three-point  sixth-order  SCCD  scheme  is  a 
promising  highly  accurate  method  for  both  derivative  computation  and  FDE  solutions. 
The  advantage  of  this  highly  accurate  three-point  staggered  scheme  is  its  potentially  wide 
application  to  geophysical  computations,  especially  to  atmospheric  and  oceanic  numerical 
models. 

(2)  This  scheme  requires  satisfaction  of  FDE  not  only  on  the  interior  grid  points,  but  also  on 
the  boundary  nodes. 

(3)  For  periodic  boundaries,  the  SCCD  scheme  has  sixth-order  accuracy  at  all  grid  points  in¬ 
cluding  boundary  nodes.  But  for  nonperiodic  boundaries,  an  additional  fifth-order  bound¬ 
ary  condition  is  needed  at  the  boundary  nodes  to  close  the  system. 

(4)  The  eigenvalue  analysis  indicates  the  stability  of  the  SCCD  scheme. 

(5)  Two  examples  (Blasius  equation  and  semispectral  primitive  equation  ocean  model)  show 
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great  reduction  of  truncation  error  using  the  SCCD  scheme.  For  the  semispectral  ocean 
model,  the  SCCD  scheme  reduces  truncation  errors  roughly  300,  36.6,  30,  8,  3  times  com¬ 
pared  to  the  second-order  (central),  ordinary  fourth-order,  compact  fourth-order,  ordinary 
sixth-order,  and  compact  sixth-order  schemes. 

(6)  Most  ocean  models  are  hydrostatically  balanced  and  on  staggered  grids.  Such  a  balance 
sets  up  a  lowest  limit,  (Ax,  A u)l,  for  the  grid  spacings.  If  the  model  resolution  is  finer 
than  that  limit,  nonhydrostatic  ocean  models  should  be  used.  If  the  hydrostatic  balance 
is  kept  unchanged,  the  grid  spacing  cannot  be  artificially  small.  A  logical  approach  is  to 
use  high-order  difference  schemes  such  as  the  SCCD  scheme. 

(7)  Future  studies  include  applying  SCCD  scheme  to  nonuniform,  staggered  grid  systems,  as 
well  as  designing  even  higher-order  schemes  such  as  eighth-order  SCCD  scheme. 


APPENDIX 

STAGGERED  FINITE  DIFFERENCE  SCHEMES 

Al.  Second-Order  Scheme 

The  staggered  second-order  scheme  is  given  by 

(dp\  _  Pi+1/2  ~  Pi-1/2  1  (d3p\  2 

\dxJi~  h  24\dx*)f  ’ 

A2.  Fourth-Order  Schemes 
Ordinary  scheme 

The  staggered  fourth-order  ordinary  scheme  is  given  by 

( &P\  Pi-3/2  ~  27pi_i/2  +  27pi_i/2 -Pi+i/2  3  /  d5p\  4 

[dxj^  24h  640  V&r5^  ’ 

9  1 

Pi  =  ^{pi-1/2  +  Pi+1/2)  -  (Pi-3/2  +  Pi+ 3/2)  +  O  ( h 4) . 

Compact  scheme 

The  staggered  fourth-order  compact  scheme  is  given  by 


Pi+1/2  -  Pi-1/2  17  / d5p\  4 

h  5, 760  \  dx5  )  { 


with 


Pi  —  Pi-1/2  +  Pi+1/2)  +  — 


h  +  O  ( h 4) . 


At  the  left  boundary,  we  use 

fdp\  (~(ll/12)p1/2+(17/24)p3/2+(3/8)p5/2  -  (5/24)p7/2  +  (l/24)p9/2)  71  f  d5p\  4 

\dxj1  h  1920  Vdx5/i  ’ 

35  35  35  7  5 

Pl  =V28Pl/2  +  32p3/2  "  64P5/2  +  32p7/2  ~  V28p9/2  +  °  ^  * 


The  right  boundary  point  has  the  similar  formulation. 
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A3.  Sixth-Order  Schemes 


Ordinary  scheme 

The  staggered  sixth-order  compact  scheme  is  given  by 


m,-- 


-9p»— 5/2  +  125^-3/2  -  2, 250pt-i/2  +  2, 250pi+i/2  -  125pi+3/2  +  9pi+5/2 


1920 h 


_  _JL_  ( h<> 

7,168  {dx5]^  ‘ 


Compact  scheme 

The  staggered  sixth-order  compact  scheme  is  given  by 


80 


with 


sfS)  +  62^)  +  9(°2.) 

\ox  J  i  \ox  J  i  \ox  ) 


_  1  — 17pi_3/2  -  189pj_i/2  +  189p<+1/2  +  17pj+3/2  61  ^  d7p' 

“  240  h  +  358, 


3,400  \dxr)i  ’ 


9  7  9 

Pi  =  ^(Pt-l/2  +  Pi+ 1/2)  +  32(Pi-3/2  +  Pi+3/2)  +  ^ 


32 


(dp\  _  ( 

\dx)  i~i  V 


dp\ 

dx)i+ 1 


A  +  0  ( h 6) . 


At  the  left  boundary,  we  use 


(dp\  ,  _  T 

V^/i  1, 


627  211  59  235  91 

Pl/2  ~  640P3/2  “  ^P5/2  +  wPV*  ~  7^9/2 


920 


48ro/*  192 

+ 


128 


443  31  3,043  fd7p\  ,6 

l,920Pll/2  960Pl3/2  +  107,520  \dx7  J  {  ’ 


Pi 


231  693  1,155  231 

T7024Pl/2  +  512P3/2  "  l7024P5/2  +  256P7/2 


495  77  21 

1  noa^9/2  ^  mo^11/2  1  no. 


-Oio/o  +0 


The  right  boundary  point  has  the  similar  formulation. 
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